Ca Notebook

Load Libraries

Load all relevant R libraries.

Connect to Relational Database

# load the PostgreSQL driver
drv <- dbDriver("Postgres")

# create a connection to the postgres database
# set the search path to the mimiciii schema
con <- dbConnect(drv, dbname = "mimic",
                 host = "localhost", port = 5432,
                 user = "mousaghannam")
dbSendQuery(con, 'set search_path to mimiciii')

Show a list of all the tables in the relational database.

# show a list of tables
dbListTables(con)

Function for more easily pulling a given table.

This is a memory conservative approaching for performing SQL queries on DB. We can pull tables in memory and do certain analyses on them, without actually having to save a given object into the R environment.

#Memory conserative approaching for performing SQL queries on DB 
#(Rewrite this function to handle multiple inputs )
tbl_mimic <- function(table) {
  table <- as.character(substitute(table))
  tbl(con, dbplyr::in_schema("mimiciii", table))
}

Create dictionary of diagnosis codes

This “dictionary” will contain all of the diagnoses events in the database associated with a given hospital admission id (hadm_id), along with a description of the icd9 code.

#Pull dictionary for diagnoses
d_icd_diagnoses <- tbl_mimic(d_icd_diagnoses) %>%
  select(icd9_code, short_title)

diagnoses_icd <-tbl_mimic(diagnoses_icd) %>%
  select(hadm_id, icd9_code, seq_num)

#ICD-9 Codes + their descriptions
icd9_dict <- diagnoses_icd %>%
  inner_join(d_icd_diagnoses, by = "icd9_code")

icd9_dict

Load Electrolyte Lab Measurements

Determine the ids corresponding to lab values for each electrolyte

Saved all of the lab_ids into R environment with collect(), for filtering in the next step.

d_labitems <- tbl_mimic(d_labitems) #Pull dictionary for lab items into memory 

d_labitems %>%
  filter(str_detect(tolower(label), "calcium")) %>%
  collect() -> lab_idsCa

lab_idsCa

Filter them

We are choosing to include the 1st and 4th row, which are strictly blood lab value draws for potassium.

#Drugs to include for potassium 
itemidsCa <- lab_idsCa$itemid[c(4)]

Pull table containing all lab measurements.

We are now taking the table labevents, which contains all the lab events in the database, and selecting the rows associated with the two blood potassium lab values itemids that we chose above.

Define repletion events within the metaview database

The d_items table in mimics contains a dictionary for all of the events associated with both the metaview and careview database. We are only going to focus on the metaview database at this point. We then want to filter this table by the label description for items that are related to potassium.

d_items <- tbl_mimic(d_items) #Pull dictionary for all lab events into memory 

d_items %>%
  filter(label %ilike% "%calcium%" ) %>%
  filter(str_detect(tolower(dbsource), "meta")) %>%
  arrange(label) -> ca_items_MV

ca_items_MV 

Although there are events sourced from either metaview or careview, these two hospital databases are actually linked to 3 different tables in the mimics database: inputevents_mv, inputevents_cv or chartevents. Since we are examining only the metavision database, it will only link to inputevents_mv and chartevents. We would like to select for only the events that link to inputevents_mv.

ca_items_MV %>% 
  filter(linksto=="inputevents_mv") %>% collect() -> ca_items_MV_IM

ca_items_MV_IM

If you look at the above table, the first row has a KCL (Bolus) with unit mL, while the 5th row has Potassium Chloride with units in mEq. According to the mimics iii website:

“Metavision records IO data using two tables: RANGESIGNALS and ORDERENTRY. These tables do not appear in MIMIC-III as they have been merged to form the > INPUTEVENTS_MV table. RANGESIGNALS contains recorded data elements which last for a fixed period of time. Furthermore, the RANGESIGNALS table recorded > information for each component of the drug separately. For example, for a norepinephrine administration there would be two components: a main order > component (norepinephrine) and a solution component (NaCl).” ¹

Therefore, we will only be selecting repletions based off of the “main order component”, and ignoring the solution. We will select for the 4th and 5th rows, which correspond to Potassiume Acetate and Potassium Chloride.

ca_items_MV_IM %>% slice(1) -> ca_items_MV_IM #ONLY KEEPING THE REPLETIONS, THE ADDITIVE AMOUNTS, NOT THE BOLUS 
ca_items_MV_IM

ca_items_MV_IM %>% pull(itemid) -> ca_items_MV_IM_vector 

Pull repletion events from metaview databse

We will now use the filtered ids for the two additive solutions above to pull electrolyte repletions from the inputevents_mv table.

Determine the amount of potassium repletions:

repEventsCa %>% count() %>% collect() -> numCaRepletions

print(paste(" There are ",pull(numCaRepletions[1,1]), " total potassium repletions that were extracted from the inputevents_mv table." ))

Join tables containing repletion events and lab events

Both the tables containing the lab draws and the repletion events have been extracted. We can now join them to begin analyzing the relationship between them.

In addition to joining the tables, we have renamed the column corresponding to the moment the fluid collection for the lab value was recorded, charttime, as charttime.lab for clarity.

repEventsCa %>% 
  inner_join(labEventsCa, by=c("subject_id" = "subject_id", "hadm_id" = "hadm_id")) %>%
  distinct() %>%
  rename(charttime.lab = charttime) %>%
  collect() -> ca_lab_repletions_MV_new

ca_lab_repletions_MV_new

Find recent repletions

We want to flag all repletion events that occurred either 24 hours before or 24 hours after a given lab value. We then filter the table, storing all of the repletions occurring BEFORE a given lab value in one table, and all of the ones occurring AFTER a given lab value in another.

#Find the RECENT Repletions
ca_lab_repletions_MV_new %>%
  mutate(charttime.lab= as_datetime(charttime.lab), endtime=as_datetime(endtime), starttime = as_datetime(starttime)) %>%
  mutate(isRecentPre = difftime(starttime, charttime.lab, units = "hours") <= 24 & difftime(starttime, charttime.lab, units = "hours") > 0 ) %>%
  mutate(isRecentPost = difftime(endtime, charttime.lab, units = "hours") >= -24 & difftime(endtime, charttime.lab, units = "hours") < 0 )  -> allRepLabEvents_ca

allRepLabEvents_ca %>% 
  filter(isRecentPre | isRecentPost)

Extract the most recent repletion, prior and after a given lab value.

We have a table that contains all of the repletion events that occur either 24 hours prior or after a given lab value. We then group the rows by each subject and hospital id, as well as the time of fluid acquisition for the lab value. Then, we create a column that generates a flag if the lab value that is chosen is the minimum starttime for all the possible starttimes. Because there are multiple lab values, as well as multiple repletion events, for a given hadm_id the table outputted has all the possible permutations. We want to select the most recent repletion AFTER a lab value, so we choose the repletion with the starttime equal to the minimum starttime for a given lab value.

#Potassium lab event PRE-repletions 
allRepLabEvents_ca %>%
  filter(isRecentPre) %>% 
  group_by(subject_id, hadm_id,charttime.lab) %>%
  mutate(isMostRecentRepletion = starttime == min(starttime)) %>%
  filter(isMostRecentRepletion) %>%
  ungroup() %>% 
  group_by(subject_id, hadm_id, starttime) %>%
  mutate(isMostRecentLabEvent = charttime.lab == max(charttime.lab)) %>%
  filter(isMostRecentLabEvent)  %>% distinct() -> pre_ca_lab_repletions_MV_new

pre_ca_lab_repletions_MV_new

For the post-repletion lab values, we are now doing the same thing, but instead finding the maximum endtime prior to a given lab value, to select the most recent repletion prior to a given lab value.

#Potassium lab event POST-repletions 
allRepLabEvents_ca %>%
  filter(isRecentPost) %>% 
  group_by(subject_id, hadm_id,charttime.lab) %>%
  mutate(isMostRecentRepletion = endtime == max(endtime)) %>%
  filter(isMostRecentRepletion) %>%
  ungroup() %>% 
  group_by(subject_id, hadm_id, endtime) %>%
  mutate(isMostRecentLabEvent = charttime.lab == min(charttime.lab)) %>%
  filter(isMostRecentLabEvent)  %>% distinct() -> post_ca_lab_repletions_MV_new

post_ca_lab_repletions_MV_new

We then combine all of the values from the “pre-repletion” and “post-repletion” tables, into one table. We create a new column preVsPost to indicate if a given event is either before or after a lab value. Conceivably, the same repletion could correspond to a

#Combine them into one dataset 
ca_pre_post_repletions_MV_new <- bind_rows(list(pre_repletion = pre_ca_lab_repletions_MV_new, post_repletion = post_ca_lab_repletions_MV_new), .id = "preVsPost")
ca_pre_post_repletions_MV_new %>% distinct()

NA

Non Repletion Events

We will now determine the set of lab values for which there was no immediate repletion afterward.

Extract non-repleted values within 24 hours

In this section, we have started by filtering all the repletions prior and after a lab value, just like for the analysis on the repletions. However, this time, if there are multiple lab values that occur before a given repletion (or after), we are taking the ones after the first one.

#Potassium lab event PRE-repletions 
allRepLabEvents_ca %>%
  filter(isRecentPre) %>% 
  group_by(subject_id, hadm_id,charttime.lab) %>%
  mutate(isMostRecentRepletion = starttime == min(starttime)) %>%
  filter(isMostRecentRepletion) %>%
  ungroup() %>% 
  group_by(subject_id, hadm_id, starttime) %>%
  mutate(isMostRecentLabEvent = charttime.lab == max(charttime.lab)) %>%
  filter(!isMostRecentLabEvent)  %>% distinct() -> nonRepletLabsPre 

#Potassium lab event POST-repletions 
allRepLabEvents_ca %>%
  filter(isRecentPost) %>% 
  group_by(subject_id, hadm_id,charttime.lab) %>%
  mutate(isMostRecentRepletion = endtime == max(endtime)) %>%
  filter(isMostRecentRepletion) %>%
  ungroup() %>% 
  group_by(subject_id, hadm_id, endtime) %>%
  mutate(isMostRecentLabEvent = charttime.lab == min(charttime.lab)) %>%
  filter(!isMostRecentLabEvent)  %>% distinct() -> nonRepletLabsPost

allNonRepletions_ca <- bind_rows(list(pre_repletion = nonRepletLabsPre, post_repletion = nonRepletLabsPost), .id = "preVsPost")

Exclusions

We now want to exclude all the rows based on certain criteria.

We have a table called hadm_id_table which contains all of the hadm_ids of patients with diagnoses that would confound our results (kidney disease, etc). (WILL ADD HOW I CREATED THIS TABLE LATER).

hadm_id_table

As a result, we could simply do a semi-join to determine how many people in each group were excluded. Note, we are doing the anti-join on the table with the most recent repletion before and after a given lab result, within a 24 hour time range. We are not doing it on the set of ALL possible repletion and or lab result events.

hadm_id_table %>%
  semi_join(ca_pre_post_repletions_MV_new, by = "hadm_id") %>% 
  select(hadm_id, .id) %>%
  group_by(.id) %>%
  summarize(n = n()) %>% 
  arrange(desc(n)) %>%
  kable()

# ca_pre_post_repletions_MV_new %>%
#   ungroup()%>%
#   semi_join(hadm_id_table, by = "hadm_id") %>% 
#   select(icustay_id, orderid, .id) %>%
#   group_by(.id) %>%
#   summarize(n = n()) %>% 
#   arrange(desc(n)) %>%
#   kable()

And here is the resultant dataset after actually excluding this time, using an anti-join instead.

ca_pre_post_repletions_MV_new %>%
  anti_join(hadm_id_table, by = "hadm_id") -> postExclusionsCa_Repleted

postExclusionsCa_Repleted

###Another Exclusions

ca_pre_post_repletions_MV_new %>%
  semi_join(hadm_id_table, by = "hadm_id")

hadm_id_table %>%
  semi_join(ca_pre_post_repletions_MV_new, by = "hadm_id")

Analysis and Visualization

We are going to do analysis here on the dataset containing repletions. We are making this distinction, as later on, we will do analysis on non-repletions.

Analysis on Repletions

Table

In this table we are just doing some summary analysis on the dataset.

postExclusionsCa_Repleted %>%
  group_by( preVsPost) %>%
  summarize(n = n(), mean = mean(valuenum, na.rm = TRUE), standard_deviation = sd(valuenum, na.rm = T)) 

`summarise()` ungrouping output (override with `.groups` argument)

Histogram

Histogram of Pre and Post Repletions

Pie Chart

Non - Repletion Analysis

allNonRepletions_ca <- bind_rows(list(pre_repletion = nonRepletLabsPre, post_repletion = nonRepletLabsPost), .id = "preVsPost")

#Run Exclusions
allNonRepletions_ca %>%
  anti_join(hadm_id_table, by = "hadm_id") -> postexclusion_nonRepletions

postexclusion_nonRepletions

Table

postexclusion_nonRepletions %>%
  group_by( preVsPost) %>%
  summarize(n = n(), mean = mean(valuenum, na.rm = TRUE), standard_deviation = sd(valuenum, na.rm = T)) 

Combining the repleted vs NonRepleted pre lab values for comparison

repletedVsNonRepleted_PreV <- bind_rows(list(repleted = postExclusionsCa_Repleted, nonRepleted = postexclusion_nonRepletions), .id = "repleteVsNon") %>% filter(preVsPost == "pre_repletion")
repletedVsNonRepleted_PreV %>% distinct()

Table

#Table
repletedVsNonRepleted_PreV %>%
  group_by(repleteVsNon,flag) %>%
  summarize(observations = n(),mean = mean(valuenum, na.rm = T), std_dev = sd(valuenum, na.rm=T)) %>%
  mutate(percentage = observations / sum(observations) * 100) 

Visualization

#Visualizing the non-repleted vs repleted lab values 
ggpar(gghistogram(data = repletedVsNonRepleted_PreV, x="valuenum", fill = "flag", add="mean",position = "identity", palette = c("#00AFBB", "#E7B800"),add_density = FALSE, bins = 80, gggtheme = theme_pubr(), xlab = "Lab Value", title = "Pre Vs Post Repletion Lab Values - Potassium", ylab = "Repletions", facet.by = "repleteVsNon")
      , xlim = c(0.5,2))

Question: How can there be histogram bars that are “abnormal” on top of bars that are white? Are they just not granular enough to differentiate them? I think there are errors in how the data was flagged. Should look more into this. #### Looking at the outliers

Repleted above the threshold.

repletedVsNonRepleted_PreV %>%
  filter(repleteVsNon == "repleted") %>%
  filter(flag == "abnormal" & valuenum > 4) %>%
  ungroup()%>%
  summarize(observations = n(),mean = mean(valuenum, na.rm = T), std_dev = sd(valuenum, na.rm=T)) 

Didn’t replete Below the threshold.

Table

postexclusion_nonRepletions %>%
  filter(flag == "abnormal" & valuenum < 4) %>%
  ungroup() %>%
  summarize(observations = n(),mean = mean(valuenum, na.rm = T), std_dev = sd(valuenum, na.rm=T)) 

Visualization of not repletiong below the threshold

postexclusion_nonRepletions %>%
  filter(valuenum < 3.5 & flag == "abnormal") %>%
  ggplot(aes(x = valuenum)) +
  geom_histogram(alpha=0.6, position="identity",bins=14, color="black") + scale_x_continuous(limits = c(1,3.5)) + theme_pubr()

###Amount of people that replete (vs Non-Replete) across lab value

repletedVsNonRepleted_PreV %>%
  group_by(valuenum, repleteVsNon) %>%
  summarize(observations = n() ) %>%
  ggplot(aes(fill=repleteVsNon, y=observations, x=valuenum)) + 
    geom_bar(position="fill", stat="identity") + scale_x_continuous()  + theme_pubr()

Time of Day Analysis

Examining the time of day for each repletion and lab value.

repletedVsNonRepleted_PreV %>%
  mutate(RepletionHour = hour(starttime), LabHour = hour(charttime.lab)) %>%
  gghistogram(x = c("LabHour","RepletionHour"), bins = 24, palette  = "uchicago", alpha = 0.5, merge = T, xlab = "Hours", ylab = "Number of Occurences")

repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab) + 4, RepletionHour = hour(starttime)) %>%
  pivot_longer(cols = c(LabHour,RepletionHour), values_to="Hour", names_to="WhichHour") -> yy

ordered(yy$WhichHour)
yy$WhichHour <- factor(yy$WhichHour, ordered = FALSE)
yy$WhichHour <- relevel(yy$WhichHour,"RepletionHour")
ggplot(yy, aes(x = Hour,fill=WhichHour)) + geom_histogram(position="identity",bins = 24,alpha=0.6) + 
  scale_x_continuous(breaks=seq(0,24,2))  + theme_pubr()

ggplot(yy, aes(x = Hour,fill=WhichHour)) +
  geom_histogram(alpha=0.4, position="identity",bins=24, color="black") + scale_x_continuous(breaks=seq(0,24,2)) + facet_wrap(vars(repleteVsNon))

Pre-Repletion vs Non-Repleted Lab Values

Below, we have the time of day for lab value draws that occur for lab values that were repleted, and those that did not directly lead to a repletion.

It seems that those leading to a repltion were more likely to occur in the morning, while those that did not directly lead to a repletion were more likely to occur in the afternoon.

Why is this occurring? Perhaps, repletions in the morning are routine, but there is no routine in performing repletions in the afternoon.

repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab)) %>%
  filter(preVsPost == "pre_repletion") %>%
  ggplot(aes(x = LabHour,fill=repleteVsNon)) + geom_histogram(position="identity",bins = 24,alpha=0.6, color = "black") + scale_x_continuous(breaks=seq(0,24,2))  + theme_pubr() + xlab("Hour of Day") + ylab("Number of Lab Values") + scale_fill_discrete(name = "", labels = c("Non-Repletion", "Repletion"))

Average pre-repletion lab value at each hour

Table

repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab)) %>%
  filter(preVsPost == "pre_repletion") %>%
  group_by(LabHour, repleteVsNon) %>%
  summarize(mean_lab_value = mean(valuenum, na.rm = T)) %>%
  pivot_wider(names_from = repleteVsNon, values_from = mean_lab_value)

Visualization

repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab)) %>%
  filter(preVsPost == "pre_repletion") %>%
  group_by(LabHour, repleteVsNon) %>%
  summarize(mean_lab_value = mean(valuenum, na.rm = T)) %>%
  ggbarplot(x="LabHour", y = "mean_lab_value", fill = "repleteVsNon", facet.by = "repleteVsNon" ) +
  scale_x_continuous(breaks=seq(0,24,2))  + theme_pubr()

repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab)) %>%
  filter(preVsPost == "pre_repletion") %>%
  group_by(LabHour, repleteVsNon) %>%
  summarize(mean_lab_value = mean(valuenum, na.rm = T)) %>%
  ggplot(aes(x=LabHour,y=mean_lab_value,group=repleteVsNon, color = repleteVsNon)) + geom_line() + geom_point() + scale_x_continuous(breaks=seq(0,24,2))  + theme_pubr()


repletedVsNonRepleted_PreV %>%
  ungroup() %>%
  mutate(LabHour = hour(charttime.lab)) %>%
  filter(preVsPost == "pre_repletion") %>%
  group_by(LabHour, repleteVsNon) %>%
  summarize(mean_lab_value = mean(valuenum, na.rm = T)) %>%
  ggplot(aes(x=LabHour,y=mean_lab_value,fill=repleteVsNon)) + geom_bar(stat="identity", position = "dodge", color = "black") + scale_x_continuous(breaks=seq(0,24,2))  + theme_pubr()

---

1. https://mimic.physionet.org/mimicdata/io/↩︎